Behavioral and ROI-based analyses

knitr::opts_chunk$set(warning=FALSE, message=FALSE)
options(scipen=0, digits=4)

Introduction

Below are code and outputs for the ROI-based analyses in the submitted manuscript: “Neural substrates for moral judgments of psychological versus physical harm”.

To see larger versions of any figure, right-click, copy image location, and paste the address to a new tab on your browser.

If you have any questions and/or comments, please email Lily Tsoi: lily [dot] tsoi [at] bc [dot] edu.

Packages

Install packages and load libraries

packages <- c("rmarkdown", "knitr", "tidyverse", "broom", "lme4", "ordinal", "lsmeans")
packages_new <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(packages_new)) install.packages(packages_new)
lapply(packages,library,character.only=T)

Data import

Data files can be found on GitHub: https://github.com/tsoices/psych-phys-harm

Analyses require the following files:

• PSYCH-PHYS_ROI_PSCs.csv
• ROI_mvpa_results.csv

Make sure these files are in the same directory.

files <- c("ROI_PSCs.csv", "ROI_MVPA.csv")
dat_names <- c("dat_psc_orig", "dat_mvpa_orig")
for(i in 1:length(files)) {
  assign(dat_names[i], read.csv(paste(params$directory, files[i], sep='/')))
}
# change the order of levels
dat_psc_orig$Group <- factor(dat_psc_orig$Group, levels=c("NT", "ASD"))
dat_psc_orig$Violation <- factor(dat_psc_orig$Violation, levels=c("PH", "PS", "N"))

Analyses are based on the following:

• Number of NT participants: 25
• Number of ASD participants: 16
• Number of total participants: 41

Behavioral results

Examining behavioral responses in the scanner

• DV: rating (1-4)
• predictors: Condition (physical, psychological, neutral), Group (NT, ASD)

Data organization

Organize behavioral data

# calculate mean rating as variable on y-axis
dat_behav <- dat_psc_orig%>%
  filter(Violation == 'PH' | Violation == 'PS' | Violation == 'N') %>%
  group_by(Subject, Violation, Group, Item, Key) %>%
  summarise(mean=mean(Key)) %>%
  droplevels.data.frame(.)
dat_behav$Item <- match(dat_behav$Item, unique(sort(dat_behav$Item))) # ordering items such that it doesn't care about purity items

Ratings by condition and group

ggplot(dat_behav, aes(x=Violation, y=mean, color=Violation)) +
  stat_summary(fun.data="mean_cl_boot", position=position_dodge(0.2), size=1) +
  ylim(1,4) +
  facet_wrap(~Group, ncol=2, labeller=labeller(Group=c(NT="Neurotypical", ASD="ASD")), scales="free_y") +
  scale_x_discrete(labels=c('Physical','Psychological', 'Neutral')) +
  scale_colour_manual(name="Condition", labels=c("Physical", "Psychological", "Neutral"), values=c("red", "darkorchid4", "slategray")) +
  ylab("Rating\n(1=not at all, 4=very)") +
  xlab("Condition") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=16,face="bold"),
        axis.title.y=element_text(margin=margin(r=18)),
        axis.title.x=element_text(margin=margin(t=18)),
        plot.title=element_text(size=18,face="bold", margin=margin(b=20), hjust=0.5),
        legend.text=element_text(size=14),
        legend.title=element_text(size=14,face="bold"),
        strip.text=element_text(size=14),
        panel.grid.major.x=element_blank(),
        panel.grid.major.y=element_blank())

Ratings by item and group

ggplot(dat_behav, aes(y=mean, x=Item, color=Violation)) +
  stat_summary(fun.data="mean_cl_boot", na.rm=TRUE) +
  ylim(1,4) +
  facet_wrap(~Group, ncol=2, labeller=labeller(Group=c(NT="Neurotypical", ASD="ASD")), scales="free_y") +
  scale_colour_manual(name="Condition", labels=c("Physical", "Psychological", "Neutral"), values=c("red", "darkorchid4", "slategray")) +
  ylab("Rating\n(1=not at all, 4=very)") +
  xlab("Item") +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=16,face="bold"),
        axis.title.y=element_text(margin=margin(r=18)),
        axis.title.x=element_text(margin=margin(t=18)),
        plot.title=element_text(size=18,face="bold", margin=margin(b=20), hjust=0.5),
        legend.text=element_text(size=14),
        legend.title=element_text(size=14,face="bold"),
        strip.text=element_text(size=14),
        panel.grid.major.x=element_blank(),
        panel.grid.major.y=element_blank())

Analyses

Define the model

dat_behav$Key <- factor(dat_behav$Key)
model_behav <- clmm(Key ~ Violation*Group + (1|Subject) + (1|Item), 
                    data=dat_behav,
                    link="probit",
                    na.action=na.omit)

Test interaction between Group and Condition

anova(model_behav, update(model_behav, . ~ . -Violation:Group))

Likelihood ratio tests of cumulative link models:
 
                                             formula:                                             link:  threshold:
update(model_behav, . ~ . - Violation:Group) Key ~ Violation + Group + (1 | Subject) + (1 | Item) probit flexible  
model_behav                                  Key ~ Violation * Group + (1 | Subject) + (1 | Item) probit flexible  

                                             no.par  AIC logLik LR.stat df Pr(>Chisq)
update(model_behav, . ~ . - Violation:Group)      9 2443  -1213                      
model_behav                                      11 2445  -1212    2.07  2       0.36

Test main effect of condition

anova(update(model_behav, . ~ . - Violation:Group), update(model_behav, . ~ . - Violation:Group - Violation))

Likelihood ratio tests of cumulative link models:
 
                                                         formula:                                             link:  threshold:
update(model_behav, . ~ . - Violation:Group - Violation) Key ~ Group + (1 | Subject) + (1 | Item)             probit flexible  
update(model_behav, . ~ . - Violation:Group)             Key ~ Violation + Group + (1 | Subject) + (1 | Item) probit flexible  

                                                         no.par  AIC logLik LR.stat df Pr(>Chisq)    
update(model_behav, . ~ . - Violation:Group - Violation)      7 2515  -1251                          
update(model_behav, . ~ . - Violation:Group)                  9 2443  -1213      76  2     <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

lsmeans(model_behav, pairwise ~ Violation)

NOTE: Results may be misleading due to involvement in interactions

$lsmeans
 Violation  lsmean     SE df asymp.LCL asymp.UCL
 PH         1.0778 0.1362 NA    0.8108    1.3448
 PS         0.7947 0.1350 NA    0.5300    1.0594
 N         -0.9417 0.1421 NA   -1.2203   -0.6631

Results are averaged over the levels of: Group 
Confidence level used: 0.95 

$contrasts
 contrast estimate     SE df z.ratio p.value
 PH - PS    0.2831 0.1347 NA   2.102  0.0895
 PH - N     2.0195 0.1534 NA  13.163  <.0001
 PS - N     1.7364 0.1519 NA  11.432  <.0001

Results are averaged over the levels of: Group 
P value adjustment: tukey method for comparing a family of 3 estimates 

ROI-based univariate results

Data organization

Analyses described here are over the entire time course. Future versions of this document will let you select a time window and automatically refresh the outputs related to that time window.

Time points of interest (seconds):

• entire time course: 6-26
• background: 6-10
• action: 12-14
• outcome: 16-18
• intent: 20-22
• judgment: 24-26

time_entire <- dat_psc_orig %>% filter(Timepoint >= 6 & Timepoint <= 26)
time_background <- dat_psc_orig %>% filter(Timepoint >= 6 & Timepoint <= 10)
time_action <- dat_psc_orig %>% filter(Timepoint >= 12 & Timepoint <= 14)
time_outcome <- dat_psc_orig %>% filter(Timepoint >= 16 & Timepoint <= 18)
time_intent <- dat_psc_orig %>% filter(Timepoint >= 20 & Timepoint <= 22)
time_judgment <- dat_psc_orig %>% filter(Timepoint >= 24 & Timepoint <= 26)
dat_psc_long <- 
  time_entire %>% 
  filter(Violation == 'PH' | Violation == 'PS' | Violation == 'N') %>%
  group_by(Subject, Violation, ROI, Group, Item, Key) %>%
  summarise(PSC=mean(PSC)) %>%
  droplevels.data.frame(.)
dat_psc_long$Item <- as.factor(dat_psc_long$Item)

NT

Subset data to NT group only and define model

data_nt <- subset(dat_psc_long, Group == "NT")
model_nt <- lmer(PSC ~ 
                   Violation*ROI +
                   (1|Subject) + (1|Item), data=data_nt, REML=FALSE)
model_nt_1 <- lmer(PSC ~ 
                     Violation*ROI +
                     (Violation|Subject) + (1|Item), data=data_nt, REML=FALSE)

Test the condition x ROI interaction

anova(model_nt, update(model_nt, . ~ . - Violation:ROI))

Data: data_nt
Models:
update(model_nt, . ~ . - Violation:ROI): PSC ~ Violation + ROI + (1 | Subject) + (1 | Item)
model_nt: PSC ~ Violation * ROI + (1 | Subject) + (1 | Item)
                                        Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_nt, . ~ . - Violation:ROI)  9 5907 5961  -2945     5889                        
model_nt                                15 5918 6009  -2944     5888  0.96      6       0.99

Test the main effect of condition

anova(update(model_nt, . ~ . - Violation:ROI), update(model_nt, . ~ . - Violation:ROI - Violation))

Data: data_nt
Models:
update(model_nt, . ~ . - Violation:ROI - Violation): PSC ~ ROI + (1 | Subject) + (1 | Item)
update(model_nt, . ~ . - Violation:ROI): PSC ~ Violation + ROI + (1 | Subject) + (1 | Item)
                                                    Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
update(model_nt, . ~ . - Violation:ROI - Violation)  7 5911 5954  -2949     5897                          
update(model_nt, . ~ . - Violation:ROI)              9 5907 5961  -2945     5889  8.42      2      0.015 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# test pairwise contrasts
lsmeans(model_nt, pairwise ~ Violation)

Loading required namespace: lmerTest
NOTE: Results may be misleading due to involvement in interactions

$lsmeans
 Violation  lsmean      SE    df lower.CL upper.CL
 PH        0.11599 0.03906 52.59  0.03764   0.1943
 PS        0.15406 0.03906 52.59  0.07571   0.2324
 N         0.02426 0.03906 52.59 -0.05409   0.1026

Results are averaged over the levels of: ROI 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

$contrasts
 contrast estimate      SE    df t.ratio p.value
 PH - PS  -0.03807 0.04356 36.89  -0.874  0.6600
 PH - N    0.09173 0.04356 36.89   2.106  0.1025
 PS - N    0.12980 0.04356 36.89   2.979  0.0138

Results are averaged over the levels of: ROI 
P value adjustment: tukey method for comparing a family of 3 estimates 

Test whether including random slope of condition improves model

anova(model_nt, model_nt_1)

Data: data_nt
Models:
model_nt: PSC ~ Violation * ROI + (1 | Subject) + (1 | Item)
model_nt_1: PSC ~ Violation * ROI + (Violation | Subject) + (1 | Item)
           Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
model_nt   15 5918 6009  -2944     5888                        
model_nt_1 20 5928 6048  -2944     5888  0.55      5       0.99

Figure: time course for each condition by ROI

dat_psc_orig$ROI <- factor(dat_psc_orig$ROI, levels=c("RTPJ", "LTPJ", "PC", "DMPFC"))
dat_psc_orig$Violation <- factor(dat_psc_orig$Violation, levels=c("PS", "PH", "N"))
dat_psc_nt <- dat_psc_orig%>% 
  filter(Group == 'NT' & (Violation == 'PH' | Violation == 'PS' | Violation == 'N')) %>%
  group_by(Subject, Violation, ROI, Timepoint) %>%
  summarise(PSC=mean(PSC))
cols <- c("PS"="darkorchid4", "PH"="red", "N"="slategray")
rois <- c(RTPJ="rTPJ", LTPJ="lTPJ", PC="precuneus", DMPFC="dmPFC")
ggplot(dat_psc_nt, aes(y=PSC, x=Timepoint, color=Violation, fill=Violation)) +
  geom_smooth(na.rm=TRUE) +
  facet_wrap(~ROI, ncol=4, labeller=labeller(ROI=rois)) +
  annotate("rect", xmin=5, xmax=27, ymin=-Inf, ymax=Inf, alpha=.1) +
  scale_x_continuous(limits=c(0,28), breaks=seq(0,28,2)) +
  ylab("Percent signal change (PSC)") +
  xlab("Timepoint (s)") +
  scale_fill_manual(name="Condition\n", labels=c("Psychological harm", "Physical harm", "Neutral act"), values=cols) +
  scale_colour_manual(name="Condition\n", labels=c("Psychological harm", "Physical harm", "Neutral act"), values=cols) +
  theme_bw() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=24,face="bold"),
        axis.title.y=element_text(margin=margin(r=20)),
        axis.title.x=element_text(margin=margin(t=20)),
        legend.text=element_text(size=20),
        legend.title=element_text(size=24,face="bold"),
        legend.key.size=unit(3, "lines"),
        strip.text=element_text(size=28, face="bold"),
        panel.grid.major.x=element_blank(),
        panel.grid.major.y=element_blank())

Subset data by ROI (NT only) and define models

dat_rtpj_nt <- subset(dat_psc_long, ROI == "RTPJ" & Group == "NT")
dat_ltpj_nt <- subset(dat_psc_long, ROI == "LTPJ" & Group == "NT")
dat_pc_nt <- subset(dat_psc_long, ROI == "PC" & Group == "NT")
dat_dmpfc_nt <- subset(dat_psc_long, ROI == "DMPFC" & Group == "NT")
model_rtpj_nt <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
                   dat=dat_rtpj_nt, REML=FALSE)
model_ltpj_nt <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
                   dat=dat_ltpj_nt, REML=FALSE)
model_pc_nt <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
                   dat=dat_pc_nt, REML=FALSE)
model_dmpfc_nt <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
                   dat=dat_dmpfc_nt, REML=FALSE)

rTPJ

# test pairwise contrasts
lsmeans(model_rtpj_nt, pairwise ~ Violation)

$lsmeans
 Violation lsmean      SE    df lower.CL upper.CL
 PH        0.2007 0.05077 42.85 0.098307   0.3031
 PS        0.2216 0.05077 42.85 0.119200   0.3240
 N         0.1084 0.05077 42.85 0.006032   0.2108

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

$contrasts
 contrast estimate      SE    df t.ratio p.value
 PH - PS  -0.02089 0.04785 35.72  -0.437  0.9005
 PH - N    0.09228 0.04785 35.72   1.929  0.1454
 PS - N    0.11317 0.04785 35.72   2.365  0.0597

P value adjustment: tukey method for comparing a family of 3 estimates 

lTPJ

# test pairwise contrasts
lsmeans(model_ltpj_nt, pairwise ~ Violation)

$lsmeans
 Violation lsmean      SE    df lower.CL upper.CL
 PH        0.2595 0.06081 39.21  0.13651   0.3824
 PS        0.2707 0.06081 39.20  0.14768   0.3936
 N         0.1401 0.06081 39.21  0.01711   0.2630

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

$contrasts
 contrast estimate      SE    df t.ratio p.value
 PH - PS  -0.01118 0.07468 35.54  -0.150  0.9877
 PH - N    0.11940 0.07468 35.55   1.599  0.2594
 PS - N    0.13058 0.07468 35.54   1.748  0.2017

P value adjustment: tukey method for comparing a family of 3 estimates 

precuneus

# test pairwise contrasts
lsmeans(model_pc_nt, pairwise ~ Violation)

$lsmeans
 Violation   lsmean      SE    df lower.CL upper.CL
 PH        -0.03654 0.04516 49.59 -0.12727  0.05419
 PS         0.03257 0.04516 49.59 -0.05817  0.12330
 N         -0.10693 0.04516 49.59 -0.19767 -0.01620

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

$contrasts
 contrast estimate      SE    df t.ratio p.value
 PH - PS   -0.0691 0.04934 35.44  -1.401  0.3514
 PH - N     0.0704 0.04934 35.44   1.427  0.3382
 PS - N     0.1395 0.04934 35.44   2.827  0.0204

P value adjustment: tukey method for comparing a family of 3 estimates 

dmPFC

# test pairwise contrasts
lsmeans(model_dmpfc_nt, pairwise ~ Violation)

$lsmeans
 Violation   lsmean      SE    df lower.CL upper.CL
 PH         0.04848 0.07077 33.34 -0.09544   0.1924
 PS         0.10088 0.07077 33.34 -0.04304   0.2448
 N         -0.03575 0.07077 33.34 -0.17967   0.1082

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

$contrasts
 contrast estimate      SE  df t.ratio p.value
 PH - PS  -0.05240 0.07065 513  -0.742  0.7387
 PH - N    0.08423 0.07065 513   1.192  0.4584
 PS - N    0.13663 0.07065 513   1.934  0.1302

P value adjustment: tukey method for comparing a family of 3 estimates 

ASD

Subset data to ASD group only and define model

data_asd <- subset(dat_psc_long, Group == "ASD")
model_asd <- lmer(PSC ~ 
                   Violation*ROI +
                   (1|Subject) + (1|Item), data=data_asd, REML=FALSE)
# model with slope of condition does not converge, so no comparison will be made between model w/ slope and model w/o slope.
# model_asd_1 <- lmer(PSC ~ 
#                      Violation*ROI +
#                      (Violation|Subject) + (1|Item), data=data_asd, REML=FALSE)

Test the condition x ROI interaction

anova(model_asd, update(model_asd, . ~ . - Violation:ROI))

Data: data_asd
Models:
update(model_asd, . ~ . - Violation:ROI): PSC ~ Violation + ROI + (1 | Subject) + (1 | Item)
model_asd: PSC ~ Violation * ROI + (1 | Subject) + (1 | Item)
                                         Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_asd, . ~ . - Violation:ROI)  9 3506 3556  -1744     3488                        
model_asd                                15 3517 3600  -1744     3487  0.85      6       0.99

Test the main effect of condition

anova(update(model_asd, . ~ . - Violation:ROI), update(model_asd, . ~ . - Violation:ROI - Violation))

Data: data_asd
Models:
update(model_asd, . ~ . - Violation:ROI - Violation): PSC ~ ROI + (1 | Subject) + (1 | Item)
update(model_asd, . ~ . - Violation:ROI): PSC ~ Violation + ROI + (1 | Subject) + (1 | Item)
                                                     Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)   
update(model_asd, . ~ . - Violation:ROI - Violation)  7 3513 3551  -1749     3499                           
update(model_asd, . ~ . - Violation:ROI)              9 3506 3556  -1744     3488  10.9      2     0.0042 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# test pairwise contrasts
lsmeans(model_asd, pairwise ~ Violation)

NOTE: Results may be misleading due to involvement in interactions

$lsmeans
 Violation  lsmean      SE    df lower.CL upper.CL
 PH        0.11668 0.03348 37.28  0.04886  0.18451
 PS        0.17087 0.03348 37.28  0.10304  0.23870
 N         0.02851 0.03348 37.28 -0.03932  0.09633

Results are averaged over the levels of: ROI 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

$contrasts
 contrast estimate      SE   df t.ratio p.value
 PH - PS  -0.05419 0.04062 39.6  -1.334  0.3852
 PH - N    0.08818 0.04062 39.6   2.171  0.0889
 PS - N    0.14237 0.04062 39.6   3.505  0.0032

Results are averaged over the levels of: ROI 
P value adjustment: tukey method for comparing a family of 3 estimates 

Figure: time course by condition for each ROI

dat_psc_orig$ROI <- factor(dat_psc_orig$ROI, levels=c("RTPJ", "LTPJ", "PC", "DMPFC"))
dat_psc_orig$Violation <- factor(dat_psc_orig$Violation, levels=c("PS", "PH", "N"))
dat_psc_asd <- dat_psc_orig%>% 
  filter(Group == 'ASD' & (Violation == 'PH' | Violation == 'PS' | Violation == 'N')) %>%
  group_by(Subject, Violation, ROI, Timepoint) %>%
  summarise(PSC=mean(PSC))
cols <- c("PS"="darkorchid4", "PH"="red", "N"="slategray")
rois <- c(RTPJ="rTPJ", LTPJ="lTPJ", PC="precuneus", DMPFC="dmPFC")
ggplot(dat_psc_asd, aes(y=PSC, x=Timepoint, color=Violation, fill=Violation)) +
  geom_smooth(na.rm=TRUE) +
  facet_wrap(~ROI, ncol=4, labeller=labeller(ROI=rois)) +
  annotate("rect", xmin=5, xmax=27, ymin=-Inf, ymax=Inf, alpha=.1) +
  scale_x_continuous(limits=c(0,28), breaks=seq(0,28,2)) +
  ylab("Percent signal change (PSC)") +
  xlab("Timepoint (s)") +
  scale_fill_manual(name="Condition\n", labels=c("Psychological harm", "Physical harm", "Neutral act"), values=cols) +
  scale_colour_manual(name="Condition\n", labels=c("Psychological harm", "Physical harm", "Neutral act"), values=cols) +
  theme_bw() +
  theme(axis.text=element_text(size=16),
        axis.title=element_text(size=24,face="bold"),
        axis.title.y=element_text(margin=margin(r=20)),
        axis.title.x=element_text(margin=margin(t=20)),
        legend.text=element_text(size=20),
        legend.title=element_text(size=24,face="bold"),
        legend.key.size=unit(3, "lines"),
        strip.text=element_text(size=28, face="bold"),
        panel.grid.major.x=element_blank(),
        panel.grid.major.y=element_blank())

Subset data by ROI (ASD only) and define models

dat_rtpj_asd <- subset(dat_psc_long, ROI == "RTPJ" & Group == "ASD")
dat_ltpj_asd <- subset(dat_psc_long, ROI == "LTPJ" & Group == "ASD")
dat_pc_asd <- subset(dat_psc_long, ROI == "PC" & Group == "ASD")
dat_dmpfc_asd <- subset(dat_psc_long, ROI == "DMPFC" & Group == "ASD")
model_rtpj_asd <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
                   dat=dat_rtpj_asd, REML=FALSE)
model_ltpj_asd <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
                   dat=dat_ltpj_asd, REML=FALSE)
model_pc_asd <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
                   dat=dat_pc_asd, REML=FALSE)
model_dmpfc_asd <- lmer(PSC ~ Violation + (1|Subject) + (1|Item),
                   dat=dat_dmpfc_asd, REML=FALSE)

rTPJ

# test pairwise contrasts
lsmeans(model_rtpj_asd, pairwise ~ Violation)

$lsmeans
 Violation  lsmean      SE    df lower.CL upper.CL
 PH        0.13685 0.06231 32.74  0.01003   0.2637
 PS        0.15947 0.06231 32.74  0.03266   0.2863
 N         0.04914 0.06231 32.74 -0.07767   0.1760

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

$contrasts
 contrast estimate     SE    df t.ratio p.value
 PH - PS  -0.02263 0.0796 36.01  -0.284  0.9565
 PH - N    0.08770 0.0796 36.01   1.102  0.5192
 PS - N    0.11033 0.0796 36.01   1.386  0.3587

P value adjustment: tukey method for comparing a family of 3 estimates 

lTPJ

# test pairwise contrasts
lsmeans(model_ltpj_asd, pairwise ~ Violation)

$lsmeans
 Violation lsmean      SE    df lower.CL upper.CL
 PH        0.2884 0.08445 26.34  0.11496   0.4619
 PS        0.3231 0.08445 26.34  0.14960   0.4966
 N         0.1849 0.08445 26.34  0.01144   0.3584

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

$contrasts
 contrast estimate      SE  df t.ratio p.value
 PH - PS  -0.03464 0.08369 420  -0.414  0.9099
 PH - N    0.10352 0.08369 420   1.237  0.4321
 PS - N    0.13816 0.08369 420   1.651  0.2257

P value adjustment: tukey method for comparing a family of 3 estimates 

precuneus

# test pairwise contrasts
lsmeans(model_pc_asd, pairwise ~ Violation)

$lsmeans
 Violation    lsmean     SE    df  lower.CL upper.CL
 PH         0.006022 0.0449 27.86 -0.085982  0.09803
 PS         0.099761 0.0449 27.86  0.007757  0.19177
 N         -0.055337 0.0449 27.86 -0.147341  0.03667

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

$contrasts
 contrast estimate      SE    df t.ratio p.value
 PH - PS  -0.09374 0.04262 34.72  -2.199  0.0854
 PH - N    0.06136 0.04262 34.72   1.440  0.3322
 PS - N    0.15510 0.04262 34.72   3.639  0.0025

P value adjustment: tukey method for comparing a family of 3 estimates 

dmPFC

# test pairwise contrasts
lsmeans(model_dmpfc_asd, pairwise ~ Violation)

$lsmeans
 Violation   lsmean      SE  df  lower.CL upper.CL
 PH         0.05501 0.06449 288 -0.071909  0.18194
 PS         0.12049 0.06449 288 -0.006437  0.24741
 N         -0.04427 0.06449 288 -0.171196  0.08265

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

$contrasts
 contrast estimate     SE  df t.ratio p.value
 PH - PS  -0.06547 0.0912 288  -0.718  0.7531
 PH - N    0.09929 0.0912 288   1.089  0.5218
 PS - N    0.16476 0.0912 288   1.807  0.1692

P value adjustment: tukey method for comparing a family of 3 estimates 

NT vs ASD

Define the models

dat_rtpj <- subset(dat_psc_long, ROI == "RTPJ")
dat_ltpj <- subset(dat_psc_long, ROI == "LTPJ")
dat_pc <- subset(dat_psc_long, ROI == "PC")
dat_dmpfc <- subset(dat_psc_long, ROI == "DMPFC")
model_rtpj <- lmer(PSC ~ Violation*Group + (1|Subject) + (1|Item),
                   dat=dat_rtpj, REML=FALSE)
model_ltpj <- lmer(PSC ~ Violation*Group + (1|Subject) + (1|Item),
                   dat=dat_ltpj, REML=FALSE)
model_pc <- lmer(PSC ~ Violation*Group + (1|Subject) + (1|Item),
                   dat=dat_pc, REML=FALSE)
model_dmpfc <- lmer(PSC ~ Violation*Group + (1|Subject) + (1|Item),
                   dat=dat_dmpfc, REML=FALSE)
model_nt_vs_asd <- lmer(PSC ~ Violation*Group*ROI + (1|Subject) + (1|Item),
                   dat=dat_psc_long, REML=FALSE)

Test interaction: Condition x Group x ROI

anova(update(model_nt_vs_asd, . ~ . -Violation:Group:ROI), update(model_nt_vs_asd, . ~ . -Violation:Group:ROI - Violation:Group))

Data: dat_psc_long
Models:
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group): PSC ~ Violation + Group + ROI + (1 | Subject) + (1 | Item) + 
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group):     Violation:ROI + Group:ROI
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI): PSC ~ Violation + Group + ROI + (1 | Subject) + (1 | Item) + 
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI):     Violation:Group + Violation:ROI + Group:ROI
                                                                       Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group) 19 9420 9544  -4691     9382                        
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI)                   21 9424 9561  -4691     9382  0.16      2       0.92

Test interaction: Condition x Group

anova(update(model_nt_vs_asd, . ~ . -Violation:Group:ROI), update(model_nt_vs_asd, . ~ . -Violation:Group:ROI - Violation:Group))

Data: dat_psc_long
Models:
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group): PSC ~ Violation + Group + ROI + (1 | Subject) + (1 | Item) + 
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group):     Violation:ROI + Group:ROI
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI): PSC ~ Violation + Group + ROI + (1 | Subject) + (1 | Item) + 
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI):     Violation:Group + Violation:ROI + Group:ROI
                                                                       Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group) 19 9420 9544  -4691     9382                        
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI)                   21 9424 9561  -4691     9382  0.16      2       0.92

Subset data by ROI

rTPJ

anova(model_rtpj, update(model_rtpj, . ~ . -Violation:Group))

Data: dat_rtpj
Models:
update(model_rtpj, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
model_rtpj: PSC ~ Violation * Group + (1 | Subject) + (1 | Item)
                                            Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_rtpj, . ~ . - Violation:Group)  7 2415 2452  -1201     2401                        
model_rtpj                                   9 2419 2466  -1201     2401     0      2          1

lTPJ

anova(model_ltpj, update(model_ltpj, . ~ . -Violation:Group))

Data: dat_ltpj
Models:
update(model_ltpj, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
model_ltpj: PSC ~ Violation * Group + (1 | Subject) + (1 | Item)
                                            Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_ltpj, . ~ . - Violation:Group)  7 3095 3131  -1541     3081                        
model_ltpj                                   9 3099 3145  -1541     3081  0.04      2       0.98

precuneus

anova(model_pc, update(model_pc, . ~ . -Violation:Group))

Data: dat_pc
Models:
update(model_pc, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
model_pc: PSC ~ Violation * Group + (1 | Subject) + (1 | Item)
                                          Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_pc, . ~ . - Violation:Group)  7 1628 1665   -807     1614                        
model_pc                                   9 1631 1679   -807     1613  0.23      2       0.89

dmPFC

anova(model_dmpfc, update(model_dmpfc, . ~ . -Violation:Group))

Data: dat_dmpfc
Models:
update(model_dmpfc, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
model_dmpfc: PSC ~ Violation * Group + (1 | Subject) + (1 | Item)
                                             Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_dmpfc, . ~ . - Violation:Group)  7 1657 1690   -822     1643                        
model_dmpfc                                   9 1661 1704   -822     1643  0.06      2       0.97

Test main effect: Condition

anova(update(model_nt_vs_asd, . ~ . -Violation:Group:ROI -Violation:Group), update(model_nt_vs_asd, . ~ . -Violation:Group:ROI -Violation:Group -Violation))

Data: dat_psc_long
Models:
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group): PSC ~ Violation + Group + ROI + (1 | Subject) + (1 | Item) + 
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group):     Violation:ROI + Group:ROI
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group - Violation): PSC ~ Group + ROI + (1 | Subject) + (1 | Item) + Violation:ROI + 
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group - Violation):     Group:ROI
                                                                                   Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group)             19 9420 9544  -4691     9382                        
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group - Violation) 19 9420 9544  -4691     9382     0      0     <2e-16
                                                                                      
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group)                
update(model_nt_vs_asd, . ~ . - Violation:Group:ROI - Violation:Group - Violation) ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Subset data by ROI

rTPJ

anova(update(model_rtpj, . ~ . -Violation:Group), update(model_rtpj, . ~ . -Violation:Group -Violation))

Data: dat_rtpj
Models:
update(model_rtpj, . ~ . - Violation:Group - Violation): PSC ~ Group + (1 | Subject) + (1 | Item)
update(model_rtpj, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
                                                        Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
update(model_rtpj, . ~ . - Violation:Group - Violation)  5 2419 2445  -1204     2409                          
update(model_rtpj, . ~ . - Violation:Group)              7 2415 2452  -1201     2401  7.38      2      0.025 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

lTPJ

anova(update(model_ltpj, . ~ . -Violation:Group), update(model_ltpj, . ~ . -Violation:Group -Violation))

Data: dat_ltpj
Models:
update(model_ltpj, . ~ . - Violation:Group - Violation): PSC ~ Group + (1 | Subject) + (1 | Item)
update(model_ltpj, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
                                                        Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
update(model_ltpj, . ~ . - Violation:Group - Violation)  5 3097 3123  -1543     3087                          
update(model_ltpj, . ~ . - Violation:Group)              7 3095 3131  -1541     3081  5.71      2      0.058 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

precuneus

anova(update(model_pc, . ~ . -Violation:Group), update(model_pc, . ~ . -Violation:Group -Violation))

Data: dat_pc
Models:
update(model_pc, . ~ . - Violation:Group - Violation): PSC ~ Group + (1 | Subject) + (1 | Item)
update(model_pc, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
                                                      Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)   
update(model_pc, . ~ . - Violation:Group - Violation)  5 1636 1663   -813     1626                           
update(model_pc, . ~ . - Violation:Group)              7 1628 1665   -807     1614  12.5      2     0.0019 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

dmPFC

anova(update(model_dmpfc, . ~ . -Violation:Group), update(model_dmpfc, . ~ . -Violation:Group -Violation))

Data: dat_dmpfc
Models:
update(model_dmpfc, . ~ . - Violation:Group - Violation): PSC ~ Group + (1 | Subject) + (1 | Item)
update(model_dmpfc, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
                                                         Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
update(model_dmpfc, . ~ . - Violation:Group - Violation)  5 1660 1684   -825     1650                          
update(model_dmpfc, . ~ . - Violation:Group)              7 1657 1690   -822     1643  6.95      2      0.031 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Test main effect: Group

anova(update(model_nt_vs_asd, . ~ . -Violation:Group), update(model_nt_vs_asd, . ~ . -Violation:Group -Group))

Data: dat_psc_long
Models:
update(model_nt_vs_asd, . ~ . - Violation:Group): PSC ~ Violation + Group + ROI + (1 | Subject) + (1 | Item) + 
update(model_nt_vs_asd, . ~ . - Violation:Group):     Violation:ROI + Group:ROI + Violation:Group:ROI
update(model_nt_vs_asd, . ~ . - Violation:Group - Group): PSC ~ Violation + ROI + (1 | Subject) + (1 | Item) + Violation:ROI + 
update(model_nt_vs_asd, . ~ . - Violation:Group - Group):     Group:ROI + Violation:Group:ROI
                                                         Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)    
update(model_nt_vs_asd, . ~ . - Violation:Group)         27 9436 9612  -4691     9382                            
update(model_nt_vs_asd, . ~ . - Violation:Group - Group) 27 9436 9612  -4691     9382     0      0     <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Subset data by ROI

rTPJ

anova(update(model_rtpj, . ~ . -Violation:Group), update(model_rtpj, . ~ . -Violation:Group -Group))

Data: dat_rtpj
Models:
update(model_rtpj, . ~ . - Violation:Group - Group): PSC ~ Violation + (1 | Subject) + (1 | Item)
update(model_rtpj, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
                                                    Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_rtpj, . ~ . - Violation:Group - Group)  6 2414 2446  -1201     2402                        
update(model_rtpj, . ~ . - Violation:Group)          7 2415 2452  -1201     2401  1.03      1       0.31

lTPJ

anova(update(model_ltpj, . ~ . -Violation:Group), update(model_ltpj, . ~ . -Violation:Group -Group))

Data: dat_ltpj
Models:
update(model_ltpj, . ~ . - Violation:Group - Group): PSC ~ Violation + (1 | Subject) + (1 | Item)
update(model_ltpj, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
                                                    Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_ltpj, . ~ . - Violation:Group - Group)  6 3093 3124  -1541     3081                        
update(model_ltpj, . ~ . - Violation:Group)          7 3095 3131  -1541     3081   0.3      1       0.58

precuneus

anova(update(model_pc, . ~ . -Violation:Group), update(model_pc, . ~ . -Violation:Group -Group))

Data: dat_pc
Models:
update(model_pc, . ~ . - Violation:Group - Group): PSC ~ Violation + (1 | Subject) + (1 | Item)
update(model_pc, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
                                                  Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_pc, . ~ . - Violation:Group - Group)  6 1627 1659   -807     1615                        
update(model_pc, . ~ . - Violation:Group)          7 1628 1665   -807     1614  1.16      1       0.28

dmPFC

anova(update(model_dmpfc, . ~ . -Violation:Group), update(model_dmpfc, . ~ . -Violation:Group -Group))

Data: dat_dmpfc
Models:
update(model_dmpfc, . ~ . - Violation:Group - Group): PSC ~ Violation + (1 | Subject) + (1 | Item)
update(model_dmpfc, . ~ . - Violation:Group): PSC ~ Violation + Group + (1 | Subject) + (1 | Item)
                                                     Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
update(model_dmpfc, . ~ . - Violation:Group - Group)  6 1655 1684   -822     1643                        
update(model_dmpfc, . ~ . - Violation:Group)          7 1657 1690   -822     1643     0      1       0.95

ROI-based MVPA results

Data organization

# transform data to long format
dat_mvpa_long <- dat_mvpa_orig %>%
  gather(cond, acc, which(regexpr("_", colnames(.)) > 0)) %>%
  separate(cond, c("roi", "type"), extra="drop")

Comparing mean accuracies to chance (50%)

mvpa_res <- matrix(nrow=0, ncol=6,
                   dimnames=list(NULL, c("group","col", "t", "df", "p.value", "estimate")))
groups <- levels(dat_mvpa_orig$group)
for (cl in 3:ncol(dat_mvpa_orig)) {
  for (grp in groups) {
    temp <- dat_mvpa_orig[dat_mvpa_orig$group == grp, cl]
    mvpa_res <- rbind(mvpa_res, c(grp, colnames(dat_mvpa_orig)[cl],
                                  unlist(t.test(temp, mu=.5, alternative="greater"))[c(1, 2, 3, 6)]))
  }
}
# tidy up results
mvpa_res <- tidy(mvpa_res)
mvpa_res[,3:6] <- lapply(mvpa_res[,3:6], function(x) as.numeric(as.character(x)))
mvpa_res

Comparing mean accuracies of experimental tests vs. permutation tests

mvpa_res_exp_vs_perm <- matrix(nrow=0, ncol=8,
                               dimnames=list(NULL, c("group","col1", "col2", "t", "df", "p.value", "estimate.of.exp","estimate.of.perm")))
groups <- levels(dat_mvpa_orig$group)
for (cl in seq(3, ncol(dat_mvpa_orig), by=2)) {
  for (grp in groups) {
    temp1 <- dat_mvpa_orig[dat_mvpa_orig$group == grp, cl]
    temp2 <- dat_mvpa_orig[dat_mvpa_orig$group == grp, cl+1]
    mvpa_res_exp_vs_perm <- rbind(mvpa_res_exp_vs_perm, c(grp, colnames(dat_mvpa_orig)[cl], colnames(dat_mvpa_orig)[cl + 1],
                                                          unlist(t.test(temp1, temp2, alternative="greater"))[c(1, 2, 3, 6, 7)]))
  }
}
# tidy up results
mvpa_res_exp_vs_perm <- tidy(mvpa_res_exp_vs_perm)
mvpa_res_exp_vs_perm[,4:8] <- lapply(mvpa_res_exp_vs_perm[,4:8], function(x) as.numeric(as.character(x)))
mvpa_res_exp_vs_perm

Comparing mean accuracies of NT group vs. ASD group

mvpa_res_nt_vs_asd <- matrix(nrow=0, ncol=6,
                             dimnames=list(NULL, c("col", "t", "df", "p.value", "estimate of NT", "estimate of ASD")))
groups <- levels(dat_mvpa_orig$group)
for (cl in 3:ncol(dat_mvpa_orig)) {
  temp1 <- dat_mvpa_orig[dat_mvpa_orig$group == 'NT', cl]
  temp2 <- dat_mvpa_orig[dat_mvpa_orig$group == 'ASD', cl]
  mvpa_res_nt_vs_asd <- rbind(mvpa_res_nt_vs_asd, c(colnames(dat_mvpa_orig)[cl],
                                                    unlist(t.test(temp1, temp2, alternative="two.sided"))[c(1, 2, 3, 6, 7)]))
}
# tidy up results
mvpa_res_nt_vs_asd <- tidy(mvpa_res_nt_vs_asd)
mvpa_res_nt_vs_asd[,2:6] <- lapply(mvpa_res_nt_vs_asd[,2:6], function(x) as.numeric(as.character(x)))
mvpa_res_nt_vs_asd

Figure: classification accuracies

dat_mvpa_plot <- dat_mvpa_long
dat_mvpa_plot$roi <- factor(dat_mvpa_plot$roi, levels=c("RTPJ", "LTPJ", "PC", "DMPFC"))
dat_mvpa_plot$group <- factor(dat_mvpa_plot$group, levels=c("NT", "ASD"))
ggplot(dat_mvpa_plot, aes(y=acc, x=type)) +
  stat_summary(fun.data="mean_cl_boot", position=position_dodge(0.4)) +
  facet_grid(group ~ roi, labeller=labeller(roi=c(RTPJ="rTPJ", LTPJ="lTPJ", PC="precuneus", DMPFC="dmPFC"))) +
  geom_hline(yintercept=.50) +
  ylab("Classification accuracy") +
  scale_x_discrete("Test type\n", labels=c("exp"="Experimental", "perm"="Permutation")) + 
  theme_bw() +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=16,face="bold"),
        axis.title.y=element_text(margin=margin(r=20)),
        axis.title.x=element_text(margin=margin(t=20)),
        strip.text=element_text(size=16, face="bold"),
        panel.grid.major.x=element_blank(),
        panel.grid.major.y=element_blank())